home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YCOMU.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  9KB  |  342 lines
  1. /********************************  YCOMU.C  **********************************/
  2. /******************* UNSTABLE MANIFOLD COMMANDS AND MENU *********************/
  3. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  4.  
  5.  
  6. #include "yinclud.h"
  7.  
  8. #define STEST(x)     (strcmp(CodeName,x) == 0)
  9. #define TEST3(x,y,z)    if( STEST(x) || STEST(y) || STEST(z) )
  10. #define TEST4(w,x,y,z)  if( STEST(w) || STEST(x) || STEST(y) || STEST(z) )
  11.  
  12.  
  13. /* displacement for UL and UR commands */
  14. #define     TINY 0.00000001
  15.  
  16. int     caseUs(CodeName)    /* unstable manifolds */
  17. char    *CodeName;
  18. {        
  19.     double  miss,
  20.             newton();
  21.     int     invert(), acceptStable(), i, invertFlag = NO;
  22.     char   *pCode = CodeName;
  23.  
  24.     TEST("ab") {        /* this command is not essential; it is not
  25.                    called by anything; this is very similar to
  26.                    InterrupyL in YINTRPT.C which is called by
  27.                    interrupt 'L' */
  28.         abcd();
  29.         return(1);
  30.     }
  31.     TEST("lineinfo") {
  32.         PRINT
  33.             " LI number of iterates of line( >=0 ); now %d  \n\n"
  34.             ,UIterates);
  35.         PRINT
  36.             " LC Sets the line horizontally,length .002, centered at y[]\n");
  37.         PRINT
  38.             "Alternatives:\n  SX To set the x coordinate of the starting point; now %lf\n"
  39.             ,line_x1);
  40.         PRINT
  41.             "  SY To set the y coordinate of the starting point; now %lf\n"
  42.             ,line_y1);
  43.         PRINT
  44.             "  FX To set the x coordinate of the end point; now %lf\n"
  45.             ,line_x2);
  46.         PRINT
  47.             "  FY To set the y coordinate of the end point; now %lf        \n"
  48.             ,line_y2);
  49.         PRINT
  50.             "  RDL to draw the line and its iterates on top of old picture\n");
  51.         PRINT
  52.             "  DL to draw the line and its iterates after clearing screen\n");
  53.         PRINT
  54.             "  RDL to draw the line and its iterates on top of old picture\n");
  55.         PRINT
  56.             "  LD:  number of Dots on Line; now %d \n"
  57.             ,linedots);
  58.         PRINT
  59.             "  TAV:  draw image of line one period in length, constant phase space=y1[] \n"
  60.             );
  61.         PRINT
  62.             "  RTAV:  does same as TAV  except it restores the old picture first \n"
  63.             );
  64.         return(1);
  65.     }
  66.     TEST("ld") {
  67.         linedots = (int) Entervalue((double) linedots, CHECKSET);
  68.         return(1);
  69.     }
  70.     TEST("ulen") {
  71.         PRINT
  72.             "To set the length of the unstable manifold to be computed,          \n");
  73.         PRINT
  74.             "(= number of diagonals long)                          \n");
  75.         ULength = Entervalue(ULength, CHECKSET);
  76.         return(1);
  77.     }
  78.     TEST("upix") {
  79.         PRINT
  80.             "To speed the plotting you can plot every couple pixels of the unstable   \n");
  81.         PRINT
  82.             "manifold and you can interpolate in between.  Enter number(default=4)     \n");
  83.         numPixels = (int) Entervalue((double) numPixels, CHECKSET);
  84.         return(1);
  85.     }
  86.     TEST4 ("sl1", "rsl1", "ul1", "rul1") {
  87.         if(level == 2)
  88.             scr_clr();
  89.         if(CodeName[0] == 'r') {
  90.             pCode++;/* pointer to first non 'r' character of
  91.                    CodeName */
  92.             boot = 1;
  93.         }
  94.         if(*pCode == 's') {/* see if it is a stable manifold */
  95.             if(acceptStable() == NO)
  96.                 return(1);
  97.             invertFlag = YES;
  98.         }
  99.         else
  100.             invertFlag = NO;
  101.         miss = newton(period, 1);
  102.                 /* 1 means quasi-newton as opposed to newton;
  103.                    MISS is the amount that the trajectory
  104.                    misses being closed by */
  105.         for(i = 0; i < 8 && miss > TINY / 100.; i++)
  106.             miss = newton(period, 1);
  107.  
  108.         store(y, y + eqnsL, dim);
  109.         y[X_coord] = y[X_coord] - TINY * (X_upper - X_lower);
  110.         LineForU(period);
  111.         LineIterator(40, period);
  112.         if(invertFlag == YES)
  113.                 /* see if it is a stable manifold; if so, we
  114.                    must re-invert */
  115.             invert();
  116.         return(1);
  117.     }
  118.     TEST4 ("sl", "rsl", "ul", "rul") {
  119.                 /* left stable of unstable manifold */
  120.         if(level == 2)
  121.             scr_clr();
  122.         if(CodeName[0] == 'r') {
  123.             pCode++;/* pointer to first non 'r' character of
  124.                    CodeName */
  125.             boot = 1;
  126.         }
  127.         if(*pCode == 's') {/* see if it is a stable manifold */
  128.             if(acceptStable() == NO)
  129.                 return(1);
  130.             invertFlag = YES;
  131.         }
  132.         else
  133.             invertFlag = NO;
  134.         miss = newton(period, 1);
  135.                 /* 1 means quasi-newton as opposed to newton;
  136.                    MISS is the amount that the trajectory
  137.                    misses being closed by */
  138.         for(i = 0; i < 8 && miss > TINY / 100.; i++)
  139.             miss = newton(period, 1);
  140.  
  141.         store(y, y + eqnsL, dim);
  142.         y[X_coord] = y[X_coord] - TINY * (X_upper - X_lower);
  143.         LineForU(2 * period);
  144.         LineIterator(40, 2 * period);
  145.         if(invertFlag == YES)
  146.                 /* see if it is a stable manifold; if so, we
  147.                    must re-invert */
  148.             invert();
  149.         return(1);
  150.     }
  151.     TEST4 ("sr", "rsr", "ur", "rur") {
  152.         if(level == 2)
  153.             scr_clr();
  154.         if(CodeName[0] == 'r') {
  155.             pCode++;/* pointer to first non 'r' character of
  156.                    CodeName */
  157.             boot = 1;
  158.         }
  159.         if(*pCode == 's') {/* see if it is a stable manifold */
  160.             if(acceptStable() == NO)
  161.                 return(1);
  162.             invertFlag = YES;
  163.         }
  164.         else
  165.             invertFlag = NO;
  166.         miss = newton(period, 1);
  167.                 /* 1 means quasi-newton as opposed to newton;
  168.                    MISS is the amount that the trajectory
  169.                    misses being closed by */
  170.         for(i = 0; i < 8 && miss > TINY / 100.; i++) {
  171.             miss = newton(period, 1);
  172.         }
  173.         store(y, y + eqnsL, dim);
  174.         y[X_coord] = y[X_coord] + TINY * (X_upper - X_lower);
  175.         LineForU(2 * period);
  176.         LineIterator(40, 2 * period);
  177.         if(invertFlag == YES)
  178.                 /* see if it is a stable manifold; if so, we
  179.                    must re-invert */
  180.             invert();
  181.         return(1);
  182.     }
  183.     TEST4 ("sr1", "rsr1", "ur1", "rur1") {
  184.         if(level == 2)
  185.             scr_clr();
  186.         if(CodeName[0] == 'r') {
  187.             pCode++;/* pointer to first non 'r' character of
  188.                    CodeName */
  189.             boot = 1;
  190.         }
  191.         if(*pCode == 's') {/* see if it is a stable manifold */
  192.             if(acceptStable() == NO)
  193.                 return(1);
  194.             invertFlag = YES;
  195.         }
  196.         else
  197.             invertFlag = NO;
  198.         miss = newton(period, 1);
  199.                 /* 1 means quasi-newton as opposed to newton;
  200.                    MISS is the amount that the trajectory
  201.                    misses being closed by */
  202.         for(i = 0; i < 8 && miss > TINY / 100.; i++)
  203.             miss = newton(period, 1);
  204.  
  205.         store(y, y + eqnsL, dim);
  206.         y[X_coord] = y[X_coord] + TINY * (X_upper - X_lower);
  207.         LineForU(period);
  208.         LineIterator(40, period);
  209.         if(invertFlag == YES)
  210.                 /* see if it is a stable manifold; if so, we
  211.                    must re-invert */
  212.             invert();
  213.         return(1);
  214.     }
  215.     TEST("sx") {
  216.         line_x1 = Entervalue(line_x1, CHECKSET);
  217.         return(1);
  218.     }
  219.     TEST("sy") {
  220.         line_y1 = Entervalue(line_y1, CHECKSET);
  221.         return(1);
  222.     }
  223.     TEST("fx") {
  224.         line_x2 = Entervalue(line_x2, CHECKSET);
  225.         return(1);
  226.     }
  227.     TEST("fy") {
  228.         line_y2 = Entervalue(line_y2, CHECKSET);
  229.         return(1);
  230.     }
  231.     TEST("li") {        /* line iterates */
  232.         UIterates = (int) Entervalue((double) UIterates, CHECKSET);
  233.         return(1);
  234.     }
  235.     return(0);
  236. }
  237.  
  238. abcd() {
  239.     double  savea[EQUATIONS],
  240.             saveb[EQUATIONS];
  241.     int     oldCycle;
  242.  
  243.     oldCycle = cycle;
  244.     cycle = 0;
  245.  
  246.     store(savea, ya, dim);    /* ya and yb are used in line drawing */
  247.     store(saveb, yb, dim);
  248.     if(ya[X_coord] == -9999.|| yb[X_coord] == -9999.)
  249.         PRINT "ya or yb is not set \n");
  250.     else {
  251.         if(printer > 1)
  252.             PRINT "Number of images now = %d\n", images);
  253.  
  254.         if(ya[0] != -9999.&& yb[0] != -9999.) {
  255.             makeLine(savea, saveb, ya, yb, images);
  256.         }
  257.         if(ya[0] != -9999.&& yb[0] != -9999.
  258.                 && yc[0] != -9999.) {
  259.             makeLine(savea, saveb, yb, yc, images);
  260.         }
  261.         if(ya[0] != -9999.&& yb[0] != -9999.
  262.                 && yc[0] != -9999.&& yd[0] != -9999.) {
  263.             makeLine(savea, saveb, yc, yd, images);
  264.             store(ya, savea, dim);
  265.             makeLine(savea, saveb, yd, ya, images);
  266.         }
  267.     }
  268.     cycle = oldCycle;
  269. }
  270.  
  271. makeLine(savea, saveb, y1val, y2, its)
  272. double *savea,
  273.        *saveb,
  274.        *y1val,
  275.        *y2;
  276. int     its;
  277. {
  278.     store(ya, y1val, dim);
  279.     store(yb, y2, dim);
  280.     LineForB(ya, yb);    /* does not do much */
  281.     connect2 (ya[X_coord], ya[Y_coord], yb[X_coord], yb[Y_coord]);
  282.     if(its > 0)
  283.         LineIterator(its, 1);
  284.     store(ya, savea, dim);
  285.     store(yb, saveb, dim);
  286. }
  287.  
  288. ManifoldMenu() {
  289.     if(level == SETPARAM)
  290.         scr_clr();    /* in desmets pcio.a */
  291.     scr_rowcol(1, 0);
  292.  
  293.     PRINT
  294.         "            MENU for STABLE and UNSTABLE MANIFOLDS     \n\n");
  295.     PRINT
  296.         " AB:   connect ya to yb.,draw IMAGES. \n");
  297.     PRINT
  298.         " UPIX: changes the number of pixels plotted per computation; now = %d    \n\n"
  299.         ,numPixels);
  300.  
  301.     PRINT
  302.         " UL: computes left side of unstable manifold at y1[]\n");
  303.     PRINT
  304.         " UR: computes right side of unstable manifold at y1[]\n");
  305.     PRINT
  306.         " RUL and RUR are like UL & UR, but Restore previous picture first.\n\n");
  307.     PRINT
  308.         " SL,SR,RSL,RSR: the corresponding STABLE MANIFOLD commands require inverse\n"
  309.         );
  310.     PRINT
  311.         "    to be available; they do work for Henon and Differential Equations\n");
  312.     PRINT
  313.         "\nRelated commands in other menus:\n");
  314.     PRINT
  315.         " Nx or Qx: You must set y1 at a periodic point before invoking stable and\n");
  316.     PRINT
  317.         "       unstable manifold commands; see N menu\n\n"
  318.         );
  319.     PRINT
  320.         " NP:   sets the period of the orbit(currently %d)\n",
  321.         period);
  322.  
  323.     PRINT
  324.         " SD should be set somewhere between 0 and 1 to speed returns to screen\n");
  325.     erase_line();
  326.     PRINT
  327.         "    Currently SD is %lf; see P menu \n\n"
  328.         ,diameters[ScrnSec]);
  329. }
  330.  
  331.  
  332. int     acceptStable() {    /* returns whatever invert() returns */
  333.     if(invert() == NO) {
  334.                 PRINT
  335.                 "Stable Manifold processes require the inverse to be available;\n");
  336.         PRINT
  337.             "Stable Manifold processes work for Henon and Differential equations\n");
  338.         return(NO);
  339.     }
  340.     return(YES);
  341. }
  342.